### REPLICATION FILES ###
### KREWSON, SCHOENHERR, AND SHIEH -- "DID YOU HEAR ABOUT CLARENCE THOMAS?" ###
### RESEARCH AND POLITICS ###

### NECESSARY FILES:
# rdd.Rdata
# cumulative.RData
# preonly.Rdata
# allprepost.RData
# 6daysprepost.RData
# 2daysprepost.RData
# GoogleTrends.RData
# demos.RData
# webdata.RData

### R version: 4.3.3.

#Figure 1
library(rdd)
load("rdd.RData")
rdd_simple <- RDestimate(total~rdd_count,data=agg,cutpoint=0)
plot(rdd_simple)
title(xlab="Days since Propublica Publication",
      ylab="Percent of Sample Seeking Information about the Court")

#Figure 2
load("cumulative.RData")
par(mfrow=c(1,2))
plot(dat$rdd,dat$prop,col="white",xlab="Days since ProPublica publication",
     ylab="Cumulative Percentage of Population Interest",type="S",ylim=c(0,.15),yaxs="i")
abline(v=0,lty=2)
lines(dat$rdd,dat$prop,type="S",lwd=2,col="steelblue")
plot(dat$rdd,dat$prop_dif,col="white",xlab="Days since ProPublica Publication",
     ylab="Change in Cumulative Percentage",type="S",ylim=c(0,0.043),yaxs="i")
abline(v=0,lty=2)
lines(dat$rdd,dat$prop_dif,type="S",lwd=2,col="steelblue")
par(mfrow=c(1,1))

#Appendix A
summary(rdd_simple)

#Appendix B
library(stargazer)
load("preonly.RData")
summary(mod1<-glm(total~gender+as.numeric(educ)+as.numeric(pid7)+nonwhite+
                    as.numeric(faminc_new)+birthyr+poly(visits,1),agg,family="binomial"))
load("allprepost.RData")
summary(mod4<-glm(dif~gender+as.numeric(educ)+as.numeric(pid7)+nonwhite+
                    as.numeric(faminc_new)+birthyr+poly(visits,1),
                  family=binomial(link="logit"),agg))
load("6daysprepost.RData")
summary(mod3<-glm(dif~gender+as.numeric(educ)+as.numeric(pid7)+nonwhite+
                    as.numeric(faminc_new)+birthyr+poly(visits,1),
                  family=binomial(link="logit"),agg))
load("2daysprepost.RData")
summary(mod2<-glm(dif~gender+as.numeric(educ)+as.numeric(pid7)+nonwhite+
                    as.numeric(faminc_new)+birthyr+poly(visits,1),
                  family=binomial(link="logit"),agg))
stargazer(mod1,mod4,mod3,mod2)

#Appendix C
load("GoogleTrends.RData")
bw <- with(comb,IKbandwidth(dates,count,cutpoint=0))
rdd_simple <- RDestimate(count~dates,data=comb,cutpoint=0,bw=bw)
summary(rdd_simple)
plot(rdd_simple)
title(xlab="Days since Propublica Publication",ylab="Total Interest")

#Appendix F
library(xtable)
load("demos.RData")
df <- data.frame(variable=names(demos)[17:29],mean=NA,sd=NA,min=NA,max=NA,length=NA)
for(i in 17:29){
  df[i-16,2] <- mean(demos[,i],na.rm=T)
  df[i-16,3] <- sd(demos[,i],na.rm=T)
  df[i-16,4] <- min(demos[,i],na.rm=T)
  df[i-16,5] <- max(demos[,i],na.rm=T)
  df[i-16,6] <- sum(!is.na(demos[,i]))
}
xtable(df)

#Appendix H
load("webdata.RData")
hist(web$start_time_utc,breaks=30,main="Histogram of web engagement start times (UTC)",xlab="")
hist(web$date,breaks=30,main="Histogram of web engagements (by day)",xlab="")
hist(web$time,main="Histogram of web engagements by hour of day",xlab="")

#Appendix I
library(viridis)
preds1 <- predict(mod1,newdata=data.frame(gender="Female",
                                          educ=2,
                                          nonwhite=0,
                                          faminc_new=3,
                                          birthyr=1990,
                                          visits=1887,
                                          pid7=1:7),
                  type="response",se.fit=T)
preds2 <- predict(mod2,newdata=data.frame(gender="Female",
                                          educ=2,
                                          nonwhite=0,
                                          faminc_new=3,
                                          birthyr=1990,
                                          visits=1877,
                                          pid7=1:7),
                  type="response",se.fit=T)
preds3 <- predict(mod3,newdata=data.frame(gender="Female",
                                          educ=2,
                                          nonwhite=0,
                                          faminc_new=3,
                                          birthyr=1990,
                                          visits=1877,
                                          pid7=1:7),
                  type="response",se.fit=T)
preds4 <- predict(mod4,newdata=data.frame(gender="Female",
                                          educ=2,
                                          nonwhite=0,
                                          faminc_new=3,
                                          birthyr=1990,
                                          visits=1877,
                                          pid7=1:7),
                  type="response",se.fit=T)
preds1$lo <- preds1$fit-preds1$se.fit*qnorm(.975)
preds1$hi <- preds1$fit+preds1$se.fit*qnorm(.975)
preds2$lo <- preds2$fit-preds2$se.fit*qnorm(.975)
preds2$hi <- preds2$fit+preds2$se.fit*qnorm(.975)
preds3$lo <- preds3$fit-preds3$se.fit*qnorm(.975)
preds3$hi <- preds3$fit+preds3$se.fit*qnorm(.975)
preds4$lo <- preds4$fit-preds4$se.fit*qnorm(.975)
preds4$hi <- preds4$fit+preds4$se.fit*qnorm(.975)
cols <- viridis(3)
plot(1:7,preds1$fit,type="l",lwd=4,ylim=c(-0.002270816,0.099396190),col=cols[1],
     xaxt="n",xlab="",ylab="Percent Seeking Information about the Court",lty=2)
lines(1:7,preds4$fit,lwd=4)
polygon(x=c(1:7,7:1),
        y=c(preds1$lo,rev(preds1$hi)),
        border=NA,col=scales::alpha(cols[1],.5))
polygon(x=c(1:7,7:1),
        y=c(preds4$lo,rev(preds4$hi)),
        border=NA,col=scales::alpha(cols[2],.5))
legend("topright",legend=c("Information seeking post-scandal",
                           "Information seeking pre-scandal"),
       lty=1:2,col=rev(cols[1:2]),lwd=2)
abline(h=0)
axis(1,at=c(1.5,6.5),labels=c("Strong Democrat","Strong Republican"),tick=F)

